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ABSTRACT 

We have performed an analysis of the diffuse gamma-ray emission with the Fermi Large Area Telescope (LAT) in 
the Milky Way halo region, searching for a signal from dark matter annihilation or decay. In the absence of a robust 
dark matter signal, constraints are presented. We consider both gamma rays produced directly in the dark matter 
annihilation/decay and produced by inverse Compton scattering of the e + /e~ produced in the annihilation/decay. 
Conservative limits are derived requiring that the dark matter signal does not exceed the observed diffuse gamma-ray 
emission. A second set of more stringent limits is derived based on modeling the foreground astrophysical diffuse 
emission using the GALPR0P code. Uncertainties in the height of the diffusive cosmic-ray halo, the distribution 
of the cosmic -ray sources in the Galaxy, the index of the injection cosmic-ray electron spectrum, and the column 
density of the interstellar gas are taken into account using a profile likelihood formalism, while the parameters 
governing the cosmic -ray propagation have been derived from fits to local cosmic-ray data. The resulting limits 
impact the range of particle masses over which dark matter thermal production in the early universe is possible, and 
challenge the interpretation of the PAMELA/Fi?rm/-LAT cosmic ray anomalies as the annihilation of dark matter. 

Key words: dark matter - Galaxy: halo - gamma rays: diffuse background - methods: statistical 
Online-only material: color figures 


1. INTRODUCTION 

The nature of dark matter (DM) and its properties are still 
unknown, despite being one of the most widely investigated 
topics in contemporary fundamental physics. However, current 
and near future experiments are probing more and more of the 
parameter space predicted for the most popular type of DM 
candidates, weakly interacting massive particles (WIMPs; for 
a review see Bergstrom 2000; Bertone et al. 2004; Bertone 
2010). In particular, high-energy gamma-ray astronomy can 
be used to search for signatures of DM in the Milky Way 
and beyond. Gamma-rays are products of hadronization and 
radiative loss processes, and are therefore unavoidably emitted 
in the annihilation and decay of WIMPs. The propagation of 
gamma rays is mostly unaffected by the interstellar medium 
(ISM) and Galactic magnetic fields (GMFs), and therefore the 
data retain information on the morphology of the emission 
region (unlike, e.g., charged cosmic rays (CRs)). The Large Area 
Telescope (LAT), on board the Fermi gamma-ray observatory 
(Atwood et al. 2009), is now providing unprecedented high- 
quality gamma-ray data. 

We focus here on DM signatures in the diffuse gamma-ray 
emission as measured by Fermi-LXT. About 90% of the LAT 
photons are of diffuse origin. The Galactic component encodes 
information on the propagation and origin of CRs, distribution 
of CR sources, the ISM, magnetic and radiation fields in our 
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Galaxy, and unresolved point sources, while its extragalactic 
component provides a signature of energetic phenomena on 
cosmological scales. Both components are expected to include 
a contribution from DM annihilation/decay: the Galactic signal 
arises from the smooth DM halo around the Galactic center and 
Galactic substructures, while the extragalactic one arises from 
the signal of DM annihilation processes throughout the universe 
integrated over all redshifts. The extragalactic component is 
analyzed elsewhere (Abdo et al. 2010a). Here we focus on 
searching for a potential signal from DM annihilation/decay 
in the halo of the Milky Way (for previous work related to this 
topic, see Zhang et al. 2009; Papucci & Strumia 2010; Cirelli 
et al. 2010; Malyshev et al. 201 1 ; Baxter & Dodelson 2011; Lin 
et al. 2010; Abbasi et al. 2011). 

Due to the bright sources present in the Galactic center and 
the bright diffuse emission along the plane, it has been argued, 
e.g., in Serpico & Zaharijas (2008), that the region of the inner 
Galaxy, extending 10° -20° away from the Galactic plane, is 
promising in terms of the signal-to-background ratio (S/N). As 
an additional advantage, the constraints on the DM signal in 
that region become less sensitive to the unknown profile of the 
DM halo. In particular the S/N of cored profiles is only a factor 
~2 weaker than for the Navarro-Frenk- White (NFW) profile 
(Navarro et al. 1996) in this region. This should be compared 
to an order of magnitude of uncertainty when one considers 
the Galactic center region. The increase in S/N away from 
the plane is further emphasized for DM models in which DM 
annihilations result in a significant fraction of leptons in the 
final state. These leptons propagate in the Galaxy and produce 
high-energy gamma rays mainly through inverse Compton (IC) 


2 


The Astrophysical Journal, 761:91 (lBpp), 2012 December 20 


Ackermann et al. 


scattering on the interstellar radiation field (ISRF). By diffusing 
away from the Galactic center region, electrons produce an 
extended gamma-ray signal which further enhances the S/N 
at higher Galactic latitudes (Borriello et al. 2009b). 

Consequently, we investigate here a large region of interest 
(ROI) covering the central part of the Galactic halo, while 
masking out the Galactic plane. We test the LAT data for a 
contribution from the DM annihilation/decay signal by a fit of 
the spectral and spatial distributions of the expected photons in 
the ROI to the LAT data. In doing so, we take into account the 
most up-to-date modeling of the diffuse signal of astrophysical 
origin (Ackermann et al. 2012b), adapting it to the problem 
in question. This paper is organized in the following way. In 
Section 2 we describe our modeling of the diffuse gamma-ray 
emission and the way in which we parameterize it. Section 3 
outlines our general approach to fitting for DM signals in the 
presence of uncertainties in the astrophysical foregrounds. In 
Section 4 we describe the DM and gas maps used in this work 
while in Section 5 we define our data set and ROI. In Section 6 
we derive DM limits without modeling of the background, 
while Section 7 contains the details of the fit procedure for 
the limits which include modeling of the background. The 
results are presented in Section 8, while Section 9 contains 
further discussions on the background model uncertainties. In 
Section 10 we summarize and conclude. 

2. MODELING OF THE HIGH-ENERGY GALACTIC 
DIFFUSE EMISSION 

The Galactic diffuse gamma-ray emission is produced by the 
interaction of energetic CR electrons (CREs) and nucleons with 
the interstellar gas and radiation field. Its main components 
are photons from the decay of neutral pions produced in 
the interaction of the CR nucleons with the interstellar gas, 
bremsstrahlung of the CRE population on the interstellar gas 
and their IC scattering off the ISRF. Efficient modeling of the 
diffuse gamma-ray emission needs an accurate description of 
both the interstellar gas and radiation targets as well as the 
distribution of CRs in the Galaxy. 

The gas in the ISM consists mostly of atomic hydrogen 
(Hi) and, to a lesser extent, molecular hydrogen (H 2 ) which is 
concentrated along the Galactic plane. Ionized (H 11 ) hydrogen 
is subdominant, although it has a larger vertical scale (see 
Ackermann et al. 2012b; Moskalenko et al. 2004 for more 
details). Helium (He) is also important, being ~25% by mass 
of the ISM, and its distribution is assumed to follow that 
of interstellar hydrogen (Ackermann et al. 2012b). Velocity 
resolved radio surveys of the 21 cmhyperfine structure transition 
of Hi and corresponding surveys of the 2.6 mm CO J(1 — » 0) 
transition (using the CO density as a proxy for the H 2 density) 
are used to build maps of the interstellar gas in different annuli 
(Ackermann et al. 2012b) effectively providing a 3D 62 model of 
the gas distribution in the Galaxy. The conversion factors Xco 
between CO line intensity and H 2 column density have been 
observed to vary throughout the Galaxy (Abdo et al. 2010d). 
Total gas column density estimated from E(B — V) visual 
reddening maps (Schlegel et al. 1998) has been shown to be 
complementary to the one estimated from H 1 and CO surveys 
combined (Grenier et al. 2005). As described in Ackermann 
et al. (2012b), we take this into account by correcting the gas 
column density for each line of sight according to the value 


63 More precisely the model is only pseudo 3D due to the near-far ambiguity 
in the inner Galaxy (Ackermann et al. 2012b). 


derived from the E(B — V) map, except in the regions of high 
extinction (see Section 4.2). 

We use a 2D+1 cylindrically symmetric model (two spatial 
dimensions and the frequency dimension) of the ISRF, computed 
based on a model of the radiation emission of stellar populations 
and further reprocessing in the Galactic dust (Moskalenko et al. 
2006). 

We use the GALPR0P code (Strong et al. 2000) v54, to calculate 
the propagation and distribution of CRs in the Galaxy. 

The code is further used to create sky maps of the expected 
gamma-ray emission from the interactions of the CRs with 
the ISM and ISRF based on the models of the gas and radi- 
ation targets described above. GALPR0P approximates the CR 
propagation by a diffusion process into a cylindrical diffusion 
zone of half-height Zh and radius /?*. CREs and nuclei are in- 
jected by a parameterized distribution of CR sources. Energy 
losses, production of secondary particles in interactions, and 
re-acceleration of CRs in the ISM are taken into account (for 
details see Strong et al. 2000). Several important parameters en- 
ter the GALPR0P modeling: the distribution of CR sources, the 
half-height of the diffusive halo Zh , the radial extent of the halo 
Rh, the nucleon and electron injection spectrum, the normaliza- 
tion of the diffusion coefficient Do, the rigidity dependence of 
the diffusion coefficient S ( D(p ) = Dq{p / pa)~ s with p 0 = 4 
GV being the reference rigidity and /3 — v /c) the Alfven speed 
Va (parameterizing the strength of re-acceleration of CRs in the 
ISM via Alfven waves), and the velocity of the Galactic winds 
perpendicular to the Galactic plane V c . Following Ackermann 
et al. (2012b) we will parameterize the nucleon injection spec- 
trum as a broken power law in rigidity with y p 1 , y p 2 being the 
index of the spectrum before and after the break, respectively, 
and pb r p being the break rigidity. Similarly, the electron injec- 
tion spectrum is parameterized as a double broken power law 
with y e 1 , y e 2 , y e , 3 . being the index of the spectrum in the three 
rigidity zones and Pb r ,e,i , Pbr,e ,2 being the two breaks. 

We use in the following the results and formalism of 
Ackermann et al. (2012b) and we briefly summarize here the 
approach and results presented there. The reader is referred to 
Ackermann et al. (2012b) for more details and thorough dis- 
cussion. In that work a grid of models is considered with four 
values of Zh, two values of R/,, and four different models of CR 
source distributions (CRSDs). The CRSDs are set to correspond 
to the incompletely determined distributions of supernova rem- 
nants (SNRs) or tracers of star formation and collapse (pulsars, 
OB stars; (Case & Bhattacharya 1998; Yusifov & Kiigiik 2004). 
For each of these 32 models, four different assumptions on the 
column density of the hydrogen gas derived from its tracers are 
made for a total of 128 models. For each of the models, a fit of 
the model prediction to the local intensity of different CR nuclei 
and the B/C ratio is performed in order to fix the parameters 
(Do, Va , y P , 1 , Y P , 2 , Pbr)- Thus, the CR fit provides the injection 
spectrum for nuclei, the Alfven speed and the relation (zh , Do) 
for different values of Zh- In a second step, they also determine 
the electron injection spectrum from a fit of the model to the 
local spectrum of CRE using the diffusion parameters obtained 
in the first fit. However, we are not going to use the result of the 
second step, since we will instead fit the electron spectrum from 
gamma-ray data. As the last step in Ackermann et al. (2012b) 
an all-sky fit to the Fermi -LAT gamma-ray data is then per- 
formed to find the remaining parameters like the Xco factors. 61 


63 The procedure is iterated a few times in order to consider the feedback of 
the renormalized Xco factors in the propagation of CREs. 


3 


The Astrophysical Journal, 761:91 (lBpp), 2012 December 20 


Ackermann et al. 


Table 1 

CR Diffusion Parameters from Ackermann et al. (2012b) Used in This Work 


Parameter 




Value 



Halo height Zh (kpc) 

2 

4 

6 

8 

10 

15 

Diffusion coefficient Do (cm 2 s -1 ) 

2.7 x 10 28 

5.3 x 10 28 

7.1 x 10 28 8.3 x 10 28 

9.4 x 10 28 

1.0 x 10 29 

Diffusion index <5 

0.33 

0.33 

0.33 

0.33 

0.33 

0.33 

Alfven velocity va (km s _1 ) 

35.0 

33.5 

31.1 

29.5 

28.6 

26.3 

Nucleon injection index (low) Yp, l 

1.86 

1.88 

1.90 

1.92 

1.94 

1.96 

Nucleon injection index (high) y p t 2 

2.39 

2.39 

2.39 

2.39 

2.39 

2.39 

Nucleon break rigidity pbr,p (GV) 

11.5 

11.5 

11.5 

11.5 

11.5 

11.5 


The ISRF normalizations in various regions of the sky are also 
left as parameters free to vary in the fit, to account for possible 
uncertainties in the ISRF itself and the CRE distribution. In our 
analysis, we will also consider the ISRF uncertainties in more 
detail in Section 9. 

Thus, we use only the results from the first step of the 
analysis in Ackermann et al. (2012b) described above, but then 
allow for more freedom in certain parameters governing the CR 
distribution and astrophysical diffuse emission and constrain 
these parameters by fitting the models to the LAT gamma- 
ray data. Compared to Ackermann et al. (2012b), the main 
difference in our analysis is the use of a free CRSD whose 
parameters will be determined from the fit to the LAT data, as 
opposed to the four fixed CRSDs explored in Ackermann et al. 
(2012b). The procedure which we employ to obtain the best-fit 
CRSD is described in detail in Section 7. Another important 
difference is that we will keep the CRE source distribution and 
proton source distribution (which we will refer to as eCRSD 
and pCRSD) separate. This is justified as we do not know a 
priori if the bulk of CR protons and electrons is injected by the 
same class of sources. We report in Table 1 the CR diffusion 
and injection parameters for different values of Zh taken from 
Ackermann et al. (2012b) which we will use in the following . 114 
Note that the diffusion parameters in principle depend on the 
CRSD, but the dependence is weak and will be neglected in the 
following. There is also slight dependence on the parameters 
used to produce the gas maps (see Section 4) and the assumed 
X co distribution, which is also weak and will be neglected as 
well. Besides the free CRSDs and the scan over different values 
of z,h we will also scan electron injection spectra by varying the 
index y e< 2 while we will fix y e \ = 1.6, y e ,3 = 4 (Ackermann 
et al. 2012b), p br<eA = 2500 MV, and p br , e p = 2.2 TV. The last 
two parameters are left free to vary in the analysis performed in 
Ackermann et al. (2012b), although the fitted values typically 
differ by less than 20% for pbr.e. 1 and less than 10% pbr.e, 3 with 
respect to the values we report above. 

As a sanity check it should be also verified a posteriori that the 
flux and spectrum of local protons obtained after the gamma-ray 
fit are consistent with the experimentally observed ones, since 
the pCRSD fit could, potentially, affect them, and this would 
make inconsistent the use of the initial diffusion parameters, 
based on the observed proton spectrum. We, indeed, found that 
the proton spectrum after the second step is fully consistent 
with the observed one. We also verified that the normalization 
of the electron flux given by the gamma-ray fit is consistent with 
the observed one. This, though, is not strictly required for the 
consistency of the approach, since the CRE spectrum is not used 
to determine the diffusion parameters. 


64 More precisely, the zn = 2, 15 cases are not reported in Ackermann et al. 
(2012b). The values used here for Zh = 2, 15 are, however, obtained in the 
same way as the other zh cases. 


We also note that the nucleon injection spectra and the 
diffusion parameters are solely determined by fitting the local 
CR density, and neglecting the effect of DM. Injection of large 
quantities of nucleons from DM near the Galactic center that 
could alter the local abundances of CR nuclei (like the proton 
spectrum and the local B/C ratio, used to perform the fit) is 
in fact strongly excluded by anti-matter measurements. The 
anti-proton fraction, for example, recently measured by the 
PAMELA experiment up to ~ 1 00 GeV (Adriani et al. 2009b) 
is about 10 -4 above 10 GeV and thus constrains the DM 
contribution to be below this value. 

A drawback of using CR data is that they are usually affected 
by large systematics (for example, in the energy scale or in 
the solar modulation correction), so that the errors on the 
inferred diffusion parameters are larger than the statistical errors 
obtained from the fit (see Trotta et al. 2011 for a recent attempt 
to take into account these systematic effects). We checked that 
even quite large variations in these parameters affect our results 
only weakly. A more detailed discussion is deferred to Section 9. 

2.1. Limitations of the Model 

The model described above represents well the gamma-ray 
sky, although various residuals (at a ~30% level; Ackermann 
et al. 2012b), both at small and large scales, remain. These 
residuals can be ascribed to various limitations of the models: 
(1) imperfections in the modeling of gas and ISRF compo- 
nents, (2) simplified assumptions in the propagation setup (e.g., 
assumption of isotropy and homogeneity of the diffusion co- 
efficient), (3) unresolved point sources, which are expected to 
contribute to the diffuse emission at a level of 10% (Strong 
2007), and have not been taken into account in the modeling in 
Ackermann et al. (2012b), and (4) missing structures like Loop 
I (Casandiian et al. 2009) or the Galactic Bubbles/Lobes (Su 
et al. 2010). 

We will deal with these limitations in various ways. Structures 
like Loop I and the Galactic Bubbles appear mainly at high 
Galactic latitudes and their effects on the fitting can be limited 
using an ROI with limited extent in Galactic latitude. In the 
following we will consider an ROI in Galactic latitude, b, of 
5° < \b\ < 15°, and Galactic longitude, /, |/| ^ 80°; see 
Section 5. As for small-scale residuals, we believe that they 
are due to imperfections in the gas maps. In order to quantify 
their effect on the constraints of the DM properties, we will 
calculate the likelihood for several different assumptions on the 
gas total column density. 

Despite the various choices described above and the large 
freedom we leave in the model (CRSDs, Zh, index of the 
electron injection spectrum), we still see residuals in our ROI 
at the ±30% level and at >3er significance (see the figures 
in Section 8). Positive residuals, in particular, appear in various 
places in the ROI, especially in connection with the low Galactic 
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latitude extension of the Lobes and Loop I. Residuals at the 
same level appear also when we perform a fit in a high-latitude 
(|i>| > 40°) control region, where DM is not expected to 
contribute significantly, and thus seem to indicate the general 
level of accuracy achievable with the present modeling of the 
diffuse emission. Note that the residuals related to the Lobes 
and Loop I do not appear in the official Fermi - LAT diffuse 
model 65 since there they are explicitly modeled through the use 
of patches. Since the residuals do not seem obviously related 
to DM, as it would be, for example, in the case of a single 
strong positive residual near the Galactic center, we thus decide 
to focus in the following on setting limits on the possible DM 
signal, rather than searching for a DM signal. Given the presence 
of the residuals like those described above, we also decide to 
quote more generous limits at the 3cr and 5<r levels, as well as 
conservative limits without assumptions over the astrophysical 
background. A dedicated search for a DM signal will be reported 
in a forthcoming paper, where data in the Galactic plane and at 
low energies (< 1 GeV) also will be exploited in order to help 
separating a real DM signal from other astrophysical processes 
which may be responsible for the above residuals. 

3. OUTLINE OF THE APPROACH TO SET DM LIMITS 

Having a parameterized model of the astrophysical 
Galactic diffuse emission as described above, we could imagine 
exploring potential DM contributions by adding the DM com- 
ponent in the fit and performing a global joint fit of DM and 
the astrophysical model parameters. In practice, however, it is 
at the moment computationally challenging to perform such a 
large global fit. Thus, some additional simplifying assumption 
must be made. As already outlined in the previous section, the 
main simplification which we will adopt (which is also used in 
Ackermann et al. 2012b) is the splitting of the fit into two sepa- 
rate steps constraining some parameters based on measurements 
of the local CR intensities while constraining other parameters 
based on Fermi- LAT data. For the first step we rely entirely on 
Ackermann et al. (2012b) and we use the six different diffusion 
models (for six different Zh ) whose parameters are summarized 
in Table 1. We then use only gamma-ray data to perform the 
fit for the CRSDs, z,h, the electron index, gas maps with dif- 
ferent column densities of the interstellar gas, and DM. Some 
additional parameters which enter the gamma-ray fit are further 
introduced in Section 7.3. 

Our aim is to constrain the DM properties and treat the pa- 
rameters of the astrophysical diffuse gamma-ray background as 
nuisance parameters. Those parameters are typically correlated 
with the assumed DM content and it is thus important to con- 
sider them since they affect directly the DM fit. It is clear, for 
example, that the CRSD should have a large influence on the fit 
of the DM component. This can be seen from Figure 1, which 
shows that the gamma-ray signal produced by DM is some- 
what degenerate with the IC signal from CR sources placed in 
the inner Galaxy. Besides small morphological differences they 
mainly differ in the energy spectrum, which, however, is quite 
model dependent in the DM case. To explore the effect of un- 
certainties in the foreground modeling we use the LAT data to 
fit the CRSD for both the nuclei and CREs, as well as the CRE 
injection spectrum which directly affects the IC component. It 
is important to stress that to constrain DM in a self-consistent 
way, we will fit the CR distributions and DM at the same time, 


65 For a description see http://fermi.gsfc.nasa.gov/ssc/data/access/ 
lat/BackgroundModels.html. 


in order to take into account the degeneracy between the two. 
Fortunately, in our ROI, above and below the inner Galaxy, and 
around few GeV in energy, the it 0 component, which has dif- 
ferent morphology and is not degenerate with the DM signal, is 
dominant over IC and bremsstrahlung by about a factor of four 
to five. So, we expect approximately the same factor in improve- 
ment in DM constraints with respect to the fits that obtain limits 
without modeling this background. We will see in Sections 6 
and 8 that this expectation approximately holds. 

With the general approach above we set DM limits using the 
profile likelihood method (outlined in Section 7.1). Besides the 
approach above, we will also quote conservative upper limits 
using the data only (i.e., without performing any modeling of 
the astrophysical background). These conservative limits are 
along the lines of the work of Papucci & Strumia (2010) and 
Cirelli et al. (2010), which use a similar approach to set DM 
limits based on the first year of Fermi - LAT data. 

4. MAPS 
4.1. DM Maps 

The template maps used in the fits to model the DM 
contribution depend on the assumed DM distribution and the 
assumed annihilation/decay channel. Numerical simulations 
of the Milky Way scale halos indicate a smooth distribution 
that contains a large number of subhalos (Diemand et al. 
2007b; Springel et al. 2008a). The gamma-ray signal from the 
subhalo population is expected to dominate in the region of 
the outer halo, while in the inner <20° region of the Galaxy, 
its contribution is expected to be subdominant (Diemand et al. 
2007a; Springel et al. 2008b; Pieri et al. 2011). In our ROI 
the subhalo contribution therefore should be mild and we 
conservatively consider only the smooth component in this 
work. We parameterize the smooth DM density p with an NFW 
spatial profile (Navarro et al. 1996) and a cored (isothermal- 
sphere) profile (Begeman et al. 1991; Bahcall & Soneira 1980): 


NFW: 


P(r) = Pc, 




(1) 


Isothermal : 


Pir) = Po 


Rl + Rl 

r 2 + R 2 ' 


(2) 


These are traditional benchmark choices, as NFW is motivated 
by /V-body simulations, while cored profiles are instead moti- 
vated by the observations of rotation curves of galaxies and are 
also found in simulations of a Milky- Way-scale halo involving 
baryons (Maccio et al. 2012). The Einasto profile (Merritt et al. 
2006; Navarro et al. 2010) is emerging as a better fit to more 
recent numerical simulations, but for brevity we do not consider 
it here. It is expected that this profile should lead to DM limits 
stronger by ~30% in our ROI, with respect to a choice of an 
NFW profile (Cirelli et al. 2010). The main uncertainty in the 
DM halo profile comes from the poorly known (and modeled) 
baryonic effects. Indeed, with our choice of NFW and isother- 
mal profiles, our aim is to roughly bracket the uncertainties 
expected from the DM profile. For the local density of DM we 
take the value of po = 0.43 GeV cm -3 (Salucci et al. 2010), 66 

66 The measurement has a typical associated error bar of ±0.1 GeV cm~ 3 and 
a possible spread up to 0.2-0. 7 GeV cm~ 3 (Salucci et al. 2010; Cirelli et al. 
2011 ). 
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Figure 1. Upper panel: spatial (left) and spectral (right) distribution of gamma rays originating from the annihilation of a 250 GeV WIMP into bb. The left figure shows 
the expected intensity at E = 10 GeV for the full sky in Galactic coordinates. An NFW profile is assumed for the DM halo and a value of (cta v) = 4 x 10 23 cm 3 s - 1 
for the DM annihilation cross section. For comparison purposes typical spectra of the astrophysical emission from tt° decay and inverse Compton (IC) scattering 
are displayed in the right figure. The map also shows the boundaries of the region used to plot the average spectra of the right panel, and which we will use for the 
analysis described in this work. Central panel: same for a 250 GeV WIMP annihilating into /r + /u _ . The contribution from IC and from Final State Radiation (FSR) is 
shown separately in the spectrum and are superimposed in the spatial distribution. Lower panel: spatial (left) and spectral (right) distribution of the IC emission of an 
astrophysical CR source population distributed uniformly in Galactocentric radius within 1 kpc from the Galactic center and with a scale height of 200 pc. 

(A color version of this figure is available in the online journal.) 



|b| < 15°, |b| >5°, |l| < 80° 



and the scale radius is assumed to be R s — 20 kpc (NFW) and 
R c = 2.8 kpc (isothermal profile). The actual choice of the DM 
density profile does not have a major effect on our limits (see 
Section 8) as we do not consider the central few degrees of 
the Galaxy (where these distributions differ the most). A choice 
of a more extended core of ~5 kpc seems possible, although 


less favored by data (Bergstrom et al. 1998; see also Catena & 
Ullio 2010; Weber & de Boer 2010; Iocco et al. 2011 for further 
discussions on the po and DM profile uncertainties). With this 
choice our limits would worsen by a factor of <2. We also set 
the distance of the solar system from the center of the Galaxy 
to the value R Q = 8.5 kpc (Kerr & Lynden-Bell 1986). 
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For the annihilation/decay spectra we consider three channels 
with distinctly different signatures: annihilation/decay into the 
bb channel, into /x + /x~, and into r + r _ . In the first case gamma 
rays are produced through hadronization and pion decay. The 
resulting spectra are similar for all channels in which DM 
annihilations /decays produce heavy quarks and gauge bosons in 
the energy range considered here and is therefore representative 
for a large set of WIMP particle physics models. The choice 
of annihilation/decay into leptonic channels, provided by the 
second and third scenarios, is motivated by the PAMELA 
positron fraction (Adriani et al. 2009a) and the Fermi- LAT 
electrons-plus-positrons (Abdo et al. 2009) measurements (for 
interpretation of these measurements in terms of a DM signal 
see, e.g., Grasso et al. 2009; Meade et al. 2010; Bergstrom 
et al. 2009). In this case, gamma rays are dominantly produced 
through radiative processes of electrons, as well as through the 
Final State Radiation (FSR). 

We produce the DM maps with a version of GALPR0P slightly 
modified to implement custom DM profiles and injection 
spectra. For the prompt photons case GALPR0P integrates along 
the line of sight the DM gamma-ray emissivity given by 
Q y (r, E) = p 2 (oav) / 2m 2 x x dNy/dE in the annihilation 
case and by Q y (r, E ) = pTo/m x x dN y /dE in the decay 
case, m x being the DM particle mass, dN y /dE the gamma- 
ray annihilation/decay spectrum, (oau) the thermally averaged 
DM annihilation cross section and To = 1/r, the DM decay 
rate (i.e., the inverse of the DM lifetime). In the cases where 
the propagation of electrons produced in DM annihilation is 
relevant to the resulting gamma-ray emission, template maps 
are produced including the IC emission from DM annihilation/ 
decay-generated CREs, which have been propagated using the 
relevant set of propagation parameters. An example of the 
resulting DM maps at 10 GeV for two DM models (annihilation 
to bb and p. + p~ channels) is shown in Figure 1, for a DM mass 
of 250 GeV. 

We calculate the DM injection spectrum of electrons and 
gamma rays by using the PPPC4DMID tool described in Cirelli 
et al. (2011). This package provides interpolating functions 
calculated from a simulation of DM annihilation/decay with 
the PYTHIA (Sjostrand et al. 2008) event generator. It has 
recently been shown that for DM candidates with masses 
above the electroweak scale, the standard predictions for the 
annihilation/decay spectra are altered when account is taken 
of the production of electroweak gauge bosons from the FSR 
(Rachel rich et al. 2009; Ciafaloni et al. 2011a, 2011b). We 
include these electroweak corrections here, following Cirelli 
et al. (201 1), even though the effect on DM limits is marginal 
for our choice of energy range. 

4.2. Gamma-Ray Emission from CR Interactions 
with Interstellar Gas 

The astrophysical diffuse emission predicted by GALPR0P 
depends on assumptions about the distribution and column 
density of the interstellar gas entering the model. A significant 
uncertainty is related to the total gas column density due to 
the presence of dark gas. The dark gas contribution is estimated 
using the Schlegel, Finkbeiner, Davis (SFD; Schlegel et al. 1998) 
E(B — V) dust reddening map as a tracer of the total integrated 
gas column density (dark + radio-visible). The intensity of the 
velocity resolved Hi maps derived from the 21 cm survey data 
is then rescaled by the ratio between the total and radio-visible 
column densities to account for dark gas. This rescaling factor 
depends on the ratio between the dust and gas column densities 


and thus on the dust to H i ratio (d2HI) and the dust to CO ratio 
(d2CO). In Ackermann et al. (2012b) the dust to Hi and CO 
gas ratios are fixed through a regression procedure of the SFD 
map to the radio-derived gas maps (see Ackermann et al. 2012b 
for more details) which yields d2HI = 0.0137 x 10“ 20 mag cm 2 
and d2CO = 0.0458 mag (K km s -1 ) -1 . Here we will instead 
explore different values of the dust to gas ratios, to consider 
the possibility that this quantity is different in our ROI with 
respect to the all-sky derived value of Ackermann et al. (2012b) 
and thus, possibly, improve the residuals, especially at small 
scale, as discussed in Section 2.1. In particular, we will use six 
values of d2HI equally spaced in the range (0.0120-0.0170) x 
10- 20 mag cm 2 . The d2CO ratio will be instead fixed to the value 
0.04 mag (K km s -1 ) -1 . This is justified in the light of the fact 
that above ±5° of Galactic latitude, where we perform the fit, 
there is very little CO and we are thus not very sensitive to this 
parameter. Nevertheless, we checked the results for different 
values of d2CO in the range 0.03-0.06 mag (K km s -1 ) -1 and 
found no appreciable change in the results. 

The E(B — V ) map is not a reliable tracer of total column 
density when the dust reddening becomes very high. Thus a cut 
in E(B — V ) needs to be employed to exclude the regions of 
high Galactic extinction. In Ackermann et al. (2012b) a value of 
E(B — V) < 5 mag is found to be adequate and we will use it 
in the following. This cut excludes from the dark gas correction 
a narrow strip of few degrees along the Galactic plane, which is 
a region outside our ROI. 

The spin temperature of the Hi, T$ (see Ackermann et al. 
2012b for the detailed definition), used to extract the Hi 
maps from the radio data is also not very well known and 
introduces further uncertainties. However, since in Ackermann 
et al. (2012b) the total column density of the ISM is estimated 
from dust, the effect of T s is typically subdominant with respect 
to the dark gas correction. We will thus not consider it and fix 
Ts to the typical value 7'y = 150 K (Ackermann et al. 2012b). 

5. DATA SELECTION AND REGION OF INTEREST 

We use 24 months of LAT data starting from 2008 August 5 to 
2010 July 31, in the energy range between 1 GeV and 100 GeV. 
However, we use energies up to 400 GeV when deriving DM 
limits with no assumption on the astrophysical background, see 
the next section. As the data above 100 GeV have poor statistics 
they have little weight when performing the fit to the gamma-ray 
data to determine the background model. Instead, for the no- 
background limits, where no fitting is performed, high-energy 
points are relevant to determine the DM limits only for the 
hard FSR spectrum. We use only events classified as gamma 
rays in the P7CLEAN event selection and the corresponding 
P7CLEAN_V6 instrument response functions. 67 The data have 
been extracted and processed with the Fermi tools as described in 
Ackermann et al. (2012b). In order to minimize the contribution 
from the very bright Earth limb, we apply a maximum zenith 
angle cut of 90°. In addition, we also limit our data set to include 
only photons with an incidence angle from the instrument 
Z-axis of <72°. The events are divided into five logarithmically 
spaced energy bins. The total number of photons in our ROI is 
-350,000. We use a HEALPix 68 (Gorski et al. 2005) nside = 
64 pixelization scheme for the spatial binning, corresponding 
to a bin size of approximately 0. 9 x 0. 9. We further mask the 
point sources from the 1FGL source catalog (Abdo et al. 2010b) 


67 http://fermi.gsfc.nasa.gov/ssc/ 

68 http://healpix.jpl.nasa.gov 


7 


The Astrophysical Journal, 761:91 (lBpp), 2012 December 20 


Ackermann et al. 


to limit the impact of bright point sources on the fit. We mask 
all pixels within 1° of the location of the point sources. With 
this choice we remove ~25% of the photons, which leaves us 
with ~270,000. The 2FGL source catalog (Nolan et al. 2012) 
contains some further weak point sources, which we do not 
mask in order not to impact dramatically the size of our ROI. 
The effect of masking 2FGL sources is further considered in 
Section 8. Since we focus on the Galactic halo we choose an 
ROI limited by ±15° in Galactic latitude and ±80° in longitude, 
to focus on the region where the S/N for DM is the highest 
(Serpico & Zaharijas 2008), and to minimize the effect of the 
high-latitude structures like Loop I or the Galactic Lobes. This 
ROI also excludes the outer Galaxy where the DM signal should 
be lower. Furthermore, we mask the region \b\ < 5° along 
the Galactic plane, in order to reduce the uncertainty due to 
the modeling of the astrophysical and DM emission profiles 
discussed above. 

6. DM LIMITS WITH NO ASSUMPTION 
ON THE ASTROPHYSICAL BACKGROUND 

To set DM limits with no assumption on the astrophysical 
background we first convolve a given DM model with the 
exposure maps and point-spread function to obtain the counts 
expected from DM annihilation. The GaRDiAn (Gamma-Ray 
Diffuse ANalysis) software package (Ackermann et al. 2009, 
2012b) is used for this processing. The expected counts are 
then compared with the observed counts in our ROI and the 
upper limit is set to the minimum DM normalization, which 
gives counts in excess of the observed ones in at least one bin. 
More precisely, we set 3 a upper limits given by the requirement 
h,dm — 3^/h,dm > (Cowan et al. 2011), where n iD M 
is the expected number of counts from DM in the bin i and 
tij is the actual observed number of counts. It should be noted 
that the formula assumes a Gaussian model for the fluctuations, 
which is a good approximation given the large bin size and 
the number of counts per bin we use in this case (see below). 
The large Poisson noise present especially at high (>10 GeV) 
energies due to the limited number of counts per pixel affects 
the limits for DM masses above 100 GeV, weakening them 
somewhat. To reduce the Poisson noise, only for the present 
case of no-background modeling we choose a larger pixel size 
so to increase the number of counts per pixel. However, a very 
large pixel size would wash out the DM signal, diluting it in 
large regions, again weakening the limits. We chose the case 
with a pixel size of about 7° x7° (nside = 8) since it gives a 
reasonable compromise between the two competing factors. 69 
In this way limits typically improve by a factor of a few with 
respect to the case nside = 64. Limits for DM masses below 
100 GeV, instead, are only very weakly affected by the choice 
of the pixel size in the range 1°— 7°. Finally, again only for 
the present case of no modeling of the background, we use 
an extended energy range up to 400 GeV. This, in practice, is 
important only for the pf pT case for masses above 100 GeV and 
when we consider FSR only (since the /t + /i FSR annihilation 
spectrum is peaked near the energy corresponding to the DM 
mass and thus can be constrained only by using higher-energy 
data). For the other cases, instead, there is always significant 
gamma-ray emission below 100 GeV, either from prompt or IC 
photons and the extended energy range does not affect the limits 
appreciably. 


69 The mask is always defined (and applied! at nside = 64. After applying 
the mask the data (and the models) are downgraded to the larger pixel size. 


The limits derived from this analysis are discussed in 
Section 8. These constraints are about a factor of five worse 
than those obtained with a modeling of the background (see the 
next section), which is in agreement with the estimate made in 
Section 3. 

7. DM LIMITS WITH MODELING OF 
ASTROPHYSICAL BACKGROUND 

We derive a second set of upper limits taking into account 
a model of the astrophysical background. As described in 
Sections 2 and 3, the approach we use is a combined fit of 
DM and of a parameterized background model and we consider 
the uncertainties in the background model parameters through 
the profile likelihood method described below. 

7.1. Profile Likelihood and Grid Scanning 

For each DM channel and mass the model that describes 
the LAT data best maximizes the likelihood function, which is 
defined as a product running over all spatial and spectral bins i, 

Lk(8 dm ) = L k (6 DM , a) = max Ff Painp a, 0 DM ), (3) 

i 

where P ik is the Poisson distribution for observing n, events in 
bin i given an expectation value that depends on the parameter 
set (0 dm. <*). 0dm is the intensity of the DM component, and 
a represents the set of parameters that enter the astrophysical 
diffuse emission model as linear pre-factors to the individual 
model components (cf. Equation (5) below), while k denotes 
the set of parameters which enter in a nonlinear way. Individual 
GALPRDP models have been calculated for a grid of values in 
the k parameter space. For each family of models with the 
same set of nonlinear parameters k the profile likelihood curve 
is defined for each 0 dm as the likelihood, which is maximal 
over the possible choices of the parameters a for fixed 0 dm 
( see Rolke et al. 2005 and references therein). The notation a 
represents the conditional maximization of the likelihood with 
respect to these parameters. The linear part of the fit is performed 
with GaRDiAn, which for each fixed value of 0 dm finds the a 
parameters that maximizes the likelihood and the value of the 
likelihood itself at the maximum 70 (for details about fitting linear 
parameters see Section 7.3). However, since building the profile 
likelihood on a grid of 0 dm values is computationally expensive, 
we use an alternative approach including 0 D m explicitly in the 
set of parameters fitted by GaRDiAn. In this case GaRDiAn also 
computes the 0 dm value that maximizes the likelihood (the best- 
fit value 0dmo) and its It error estimated from the curvature of 
the log L k around the minimum. We then approximate the profile 
likelihood as a Gaussian in 0dm with mean 0dmo and width 0^0 ■ 
We have verified that this approximation works extremely well 
for a subset of cases for which we also explicitly computed the 
profile likelihood, tabulating it on a grid of 0 dm values. We will 
thus use this approximation throughout the rest of the analysis. 

In this way we end up with a set of k profiles of likelihood 
L k (9 dm), one for each combination of the nonlinear parameters. 
The envelope of these curves then approximates the final profile 
likelihood curve, L(0dm), where all the parameters linear and 


70 Technically, instead of maximizing the likelihood, GaRDiAn minimizes the 
(negative of) log-likelihood, — log L, using an external minimizer. For our 
analyses we used GaRDiAn with the Minuit (James & Roos 1975) minimizer. 
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Figure 2. Example of profile likelihood curves for four different DM annihilation/decay scenarios. Each curve refers to a particular model of the background. The 
envelope of the various curves approximates the global profile likelihood marginalized over the astrophysical uncertainties accounted for in our fitting procedure. 
The curve corresponding to the model setting the global minimum, ymin> is highlighted in red. The y scale is arbitrarily re-shifted so that the minimum value is zero. 
The green curve corresponds to the model setting the 3a upper limit (i.e., the model which is both part of the envelope profile likelihood and intersects the horizontal 
line located at +9). The upper limit is then effectively given by the x coordinate of the intersection point. The blue curve is similar, but for the 5a case (and intersects 
the horizontal line located at +25). For these three models the corresponding values of Zh , Ye, 2 > and d2HI are given in the caption. Panel description: 10 GeV DM 
particle decaying (DEC) into bb and NFW profile (upper left), 91 GeV DM particle annihilating (AN) into bb and NFW profile (upper right), 5 GeV DM particle 
decaying into r + r _ and NFW profile (lower left) and 750 GeV DM particle annihilating into r + r~, and NFW profile (lower right). 

(A color version of this figure is available in the online journal.) 


nonlinear have been included in the profile. 71 Examples of 
such final profile likelihood curves for specific DM models 
can be seen in Figure 2 and will be discussed more in detail 
in Section 7.3. 

Limits are calculated from the profile likelihood function by 
finding the 0 D M,iim values for which F(0 D M.iim)/E(0DM,max) is 
exp(— 9/2) and exp(— 25/2), for 3er and 5cr C.L. limits, re- 
spectively. This approximation is exact for Gaussian likelihood 
functions in one parameter and, due to invariance of the like- 
lihood function under reparameterization, it is most often also 
applicable to the non-Gaussian case (James 2006). For the case 
of handling nuisance parameters, this is not true a priori, but has 
been shown to give satisfactory properties for a variety of nui- 
sance parameter configurations (e.g., in Rolke et al. 2005; Abdo 
et al. 2010c; Ackermann et al. 2011). In particular, see also the 
recent search for the Higgs boson at the Large Hadron Collider, 
where 0(1 00) nuisance parameters need to be taken into account 
(ATLAS Collaboration et al. 2012). We therefore are confident 
that this approach gives the desired statistical properties, i.e., 
good coverage and discovery power, also in our analysis. 


71 We will sometime use in the following the term marginalizing although, 
typically, the term applies only within the framework of Bayesian analyses. In 
our frequentist approach it is called profiling. 


7.2. Free CR Source Distribution and Constrained Setup Limits 

In this section we introduce the first set of linear parameter, 
i.e., the coefficients defining the CRSDs. The remaining linear 
parameters will be introduced in the next section. 

As noted in Section 2, CRSDs (for example, the ones 
considered in Ackermann et al. 2012b) can be modeled from 
the direct observation of tracers of SNR, and so can be 
observationally biased. The uncertainty in the distribution of 
the tracers in the inner Galaxy is therefore large and should 
be taken into account in the derivation of the DM limits. We 
therefore fit the CRSD from the gamma-ray data, as described 
below. 

Due to the linearity of the propagation equation it is possible 
to combine solutions obtained from different CRSDs. To exploit 
this feature we define a parametric CRSD as the sum of step 
functions in Galactocentric radius R, with each step spanning a 
disjoint range in R: 

e, pCRSD(R) = c-' P 0(R - Ri)6(R i+1 - R). (4) 

i 

We choose seven steps with boundaries: R , = 0, 1.0, 3.0, 5.0, 
7.0, 9.0, 12.0, 20.0 kpc. The expected gamma-ray all-sky emis- 
sion for each of the 14 single-step primary e and p distributions 
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are calculated with GALPROP. It is also worth noting that a dif- 
ferent GALPRQP run needs to be done for each set of values of the 
nonlinear parameters, since, for a given e. pCRSD the output 
depends on the entire propagation setup. For more accurate out- 
put, especially in the inner Galaxy, which we are interested in, 
GALPROP is run with a finer grid in Galactocentric radius R with 
dr = 0.1 kpc, compared to the standard grid of dr = 1 kpc. 
The coefficients c e f p are set to unity for the individual GALPROP 
runs and then fitted from the gamma-ray data as described 
below. 

In order to have conservative and robust limits we constrain 
the parameter space defined above by setting c x p = c{ p — 0, 
i.e., setting to zero the e, pCRSDs in the inner Galaxy region, 
within 3 kpc of the Galactic center. In this way, potential 
e and p CR sources which would be required in the inner 
Galaxy will be potentially compensated by DM, producing 
conservative constraints. A second important reason to set the 
inner e, pCRSDs to zero is the fact that they are strongly 
degenerate with DM (especially the inner eCRSD, see Figure 1 ). 
Besides slight morphological differences, an astrophysical CRE 
source in the inner Galaxy is hardly distinguishable from a DM 
source, apart, perhaps, from differences in the energy spectrum. 
To break this degeneracy we would need to use data along the 
Galactic plane (within ±5° in latitude) since these are expected 
to be the most constraining for the e, pCRSDs in the inner 
Galaxy. However, the Galactic center region is quite complex 
and modeling it is beyond the scope of the current paper. We 
therefore defer such a study to follow-up publications. 

7.3. Fitting Procedure 

In the fit of the expected gamma-ray emission to the Fermi- 
LAT data, we determine the normalizations of the contributions 
from DM and from each step of the CRSD function defined 
in Equation (4) that best fit the data. To achieve this, we need 
to split each contribution into several components correspond- 
ing to the type of target and physical process responsible for 
the emission. The emission from n° decay depends only on 
the distribution of the CR nuclei sources, while the emission 
from bremsstrahlung and IC depends only on the distribution 
of the CRE sources (emission from interactions of secondary 
electrons produced in CR nuclei interactions is negligible above 
1 GeV). The gamma-ray emission arising from interactions of 
CRs with molecular gas traced by CO depends further on the 
assumed conversion factor X co between the CO intensity and 
the column density of the molecular gas. This conversion factor 
is uncertain and we vary it freely for each annulus. We deter- 
mine effective X co factors implicitly in the fit by splitting the 
calculated expected gamma-ray emission from CR interactions 
with molecular gas into Galactocentric annuli, which are sepa- 
rately normalized. Additionally, an isotropic component arising 
from the extragalactic gamma-ray background and misclassified 
charged particles needs to be included to fit the Fermi-L AT data. 
We do not include sources in the fit as we use a mask to filter 
the 1FGL point sources (cf. Section 5). To rule out the possi- 
bility that some bright sources might leak out of the mask and 
bias the fit, we performed test fits including explicitly the 1FGL 
point sources as a further template map, finding that the inclu- 
sion of the point sources introduces only a negligible change 
in the results. Equation (5) summarizes how we parameterize 
the expected gamma-ray emission I in the fit based on the com- 
ponents mentioned above. Each component is calculated using 
GALPROP and is available as a template map after the GALPROP 


run. In summary, the various GALPROP outputs are combined as 

< = E,{<f (ff'o +£;4 o<o) + Cf (Xemss +£,*CO ^iremss + ^ ) } 
+ «x(Z] / + Xic) + Em aiGB.mlGB'". (5) 

The sum over i is the sum over all step-like CRSD functions, 
the sum over j corresponds to the sum over all Galactocentric 
annuli (details of the procedure of a placement of the gas 
in Galactocentric annuli and their boundaries are given in 
Ackermann et al. 2012b). // denotes the gamma-ray emission 
from atomic and ionized interstellar gas while H 2 is the one 
from molecular hydrogen and IC is the IC emission. Xy and 
Xic are the prompt and IC (when present) DM contribution and 
a x the overall DM normalization. IGB m denote the Isotropic 
Gamma-ray Background (IGB) intensity for each of the five 
energy bins over which the index m runs. For better stability of 
the fit the template for IGB'" is built starting from an IGB with 
a power-law spectrum and normalization as given in Abdo et al. 
(2010e). In this way the fit coefficients aiGB.m are typically of 
the order of 1 . In all of the rest of the expression in Equation (5) 
the energy index m is implicit since we do not allow for the 
freedom of varying the GALPROP output from energy bin to 
energy bin. Finally, it should be also noted that in our case, 
where we mask ±5° along the plane, the above expression 
actually simplifies considerably since only the local ring Xco 
factor enters the sum, since all the other H 2 rings do not extend 
further than 5°from the plane. Also to be noted is the fact that, 
since in Equation (5) H 2 denotes a gamma-ray emission map, 
the expression has been already intrinsically multiplied by an 
Xco factor to convert the CO line intensity into an H 2 column 
density. We in fact normalize all the H 2 gamma-ray maps using 
the value Xco = 1 x 10 20 cm -2 (K km s -1 ) -1 . The Xco values 
in Equation (5) are thus adimensional ratios with respect to the 
reference value 1 x 1 0 20 cm~ 2 (K km s” 1 )“ 1 . With a slight abuse 
of notation we denote them also as Xco factors. 

The above expression predicts the expected gamma-ray 
counts in terms of the parameters (c‘ i ' p , X J co , anda x fora 

total of 7+7+1+5+1 = 21 parameters). GaRDiAnis used to build 
the profile likelihood for the intensity of the DM component 
a x by finding the set of parameter values that maximize the 
likelihood for a given a x . 

The outlined procedure is then repeated for each set of 
values of the nonlinear propagation and injection parameters 
to obtain the full set of profile likelihood curves. We scan over 
the following three parameters: the half-height of the diffusive 
zone Zh, the index of the electron injection spectrum y e 2 , and 
the dust-to-H 1 ratio d2HI. Specifically, we choose six values of 
Zh = 2, 4, 6, 8, 10, 15 kpc, eight values of y e . 2 linearly spaced 
between 1.925 and 2.8, and six values of d2HI linearly spaced 
in the range (0.0120-0.0170) xlO -20 mag cm 2 . Taking into 
account the seven step functions used for the e, pCRSDs we 
scan over a grid of 7x5x8 x 6 = 1680 GALPROP models 
(or rather GALPROP runs since combinations of the steps are 
effectively a single GALPROP model). 

In order to follow more easily the entire fitting procedure we 
report in Table 2 a summary of all the parameters employed in 
our analysis, linear and nonlinear, together with their range of 
variation in the fit or discrete values used in the grid. 

Figure 2 shows some examples of the profile likelihoods for 
selected DM masses and annihilation channels. The limits are 
set by first finding the absolute minimum and then looking at 
the intersection between the envelope of the various parabolae 
and the 3a and 5a horizontal lines. An important point to 
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0.012 0.013 0.014 0.015 0.016 0.017* 

d2HI [ 1 O' 20 mag cm 2 l 



Figure 3. Profile likelihood curves for zh , Ye,!, and d2Hi. The various curves refer to the case of no DM or different DM models (see the legend in the figure, where 
we mark a dominant decay (DEC) or annihilation (AN) channel and the assumed DM profile). All minima are normalized to the same level. Horizontal dotted lines 
indicate, as in Figure 2, a difference in — 2AlogL from the minimum of 9 (3cr) and 25 (5<x). 

(A color version of this figure is available in the online journal.) 


Table 2 

Summary Table of the Parameters Varied in the Fit 


Nonlinear Parameters 

Symbol 

Grid Values 

Index of the injection CRE spectrum 

Ye , 2 

1.925, 2.050, 2.175, 2.300, 2.425, 2.550, 2.675, 2.800 

Half height of the diffusive halo a 

Zh 

2, 4, 6, 8, 10, 15 kpc 

Dust to H i ratio 

d2HI 

(0.012, 0.013, 0.014, 0.015, 0.016, 0.017) x 10^ 20 mag cm 2 

Linear Parameters 

Symbol 

Range of Variation 

eCRSD and pCRSD coefficients 

e P 
c . c ■ 

0,+oo 

Local H 2 to CO factor 

v loc 
A CO 

(0-50) x 10 20 cm' 2 (K km s” 1 )-* 

IGB normalization in various energy bins 

GfiGB.m 

Free 

DM normalization 

“x 

Free 


Notes. The top part of the table shows the nonlinear parameters and the grid values at which the likelihood is computed. The bottom part shows the 
linear parameters and the range of variation allowed in the fit. The coefficients of the CRSDs are forced to be positive, except for c‘ ,p and Cj v . which 
are set to zero. The local Xco ratio is restricted to vary in the range (0-50) x 10 20 cm -2 (K km s -1 ) -1 , while cciGB.m and a x are left free to assume 
both positive and negative values. See the text for more details. 

a The parameters Do, S, va, Yp, 1 - Yp, i > Pbr,p are varied together with Zh as indicated in Table 1. 


note is that, for each DM model, the global minimum that we 
found lies within the 3(5)er regions of many different models. 
This is a basic sanity check against a bias in our procedure, 
as would be suspected if the model giving the minimum was 
inconsistent with the bulk of the other models considered. 
This point is further illustrated in Figure 3, where the profile 
likelihoods for the three nonlinear parameters, Zh, Ye, 2 , and d2HI, 
are shown. To ease the readability of the figure the profiling 


is actually performed by further grouping DM models with 
different DM masses, but keeping the different DM channels, 
DM profiles and the annihilation/decay cases separately. The 
curve for the fit without DM is also shown for comparison. 
Each resulting curve has been further rescaled to a common 
minimum, since we are interested in showing that several 
models are within — 2AlogL ^ 25 around the minimum for 
each DM fit. The Ye ,2 profile, for example, indicates that all 
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models with y e i from 1.9 to 2.4 are within — 2AlogL ^ 25 
around the minimum illustrating that the sampling around each 
of the minima for the six DM models is dense. Similarly, the 
d2HI profile indicates that all models with d2HI in the range 
(0. 1 20-0. 1 60) x 1 0 20 mag cm 2 are within 5 a from the minima 
for each of the six DM models. Finally, the Zh profile indicates 
that basically all the considered values of Zh are close to the 
absolute minima. This last result is not surprising since, within 
our low-latitude ROI, we have little sensitivity to different Zh 
values and basically all of them fit equally well. There is some 
tendency to favor higher values of zi, when DM is not included in 
the fit, while with DM the trend is inverted; although the feature 
is not extremely significant it is potentially very interesting. 

As explained in Section 3, in our analysis the DM parameter 
is the one of primary interest and we thus treat the parameters 
(linear and nonlinear) of the diffuse emission as nuisance 
parameters, which we include to take into account degeneracies 
with DM (i.e., to marginalize over them) and establish more 
robust limits. This is a reasonable assumption since, for example, 
most of the IC emission and thus the sensitivity to y e 2 comes 
from data within 5° of the Galactic plane, which are not 
considered here. Similarly, to be sensitive to Zh, high Galactic 
latitude data should be included. These plots thus should be 
regarded only as indicative of the achievable constraints, while 
a careful analysis will be deferred to later publications. 

An issue that we are not addressing here is whether or not DM 
is required in the fit, and, in the former case, finding the best 
model among different DM models. Since we have seen that 
systematic uncertainties related to the limitations in modeling 
astrophysical contributions to the Galactic diffuse emission are 
comparable in size to any DM signal we fit, we have focused 
on setting constraints on potential DM contributions to the 
Galactic diffuse emission. The systematic uncertainties in the 
Galactic diffuse emission modeling likely could be reduced by 
including the Galactic plane/Galactic center data in the analysis. 
Furthermore, a realistic study of the problem would require also 
considering other possible components that might be present 
in the halo, like contributions from a population of unresolved 
pulsars or the emission from the bubbles/lobes. We defer this 
analysis to a subsequent study. 

8. RESULTS 

Despite the various conservative choices described above, the 
resulting limits are quite stringent. Upper limits on the velocity 
averaged annihilation cross section into various channels are 
shown in Figure 4, for an NFW and isothermal profile of the DM 
halo. The limits obtained without modeling of the astrophysical 
background are compatible with the result of similar analyses 
presented in Papucci & Strumia (2010) and Cirelli et al. (2010). 
Limits with model of the background, instead, are significantly 
improved with respect to the above ones. They are competitive 
with respect to the limits from LAT searches for a signal from 
DM annihilation/decay in dwarf galaxies (Ackermann et al. 
201 1), and Galaxy clusters (Ackermann et al. 2010). 

In particular, as shown in Figure 4 for masses around 20 GeV 
the thermal relic value of the annihilation cross section is 
reached, both for the bb and r + r” channels. The limits are 
also improved over the ones derived in the analysis of dwarf 
galaxies (Ackermann et al. 2011), which did not consider the IC 
emission (in dwarfs this component is quite uncertain) and also 
improved over constraints imposed on DM annihilations from 
the absence of a measurable effect on CMB anisotropies (Galli 
et al. 2011). 


A limitation of our constraints is their dependence on poorly 
determined properties of the Milky Way DM halo, in particular 
on po, from which the normalization of the DM signal, and thus 
the limits, depends quadratically in the annihilation case and 
linearly in the decay case. We use the recent determination 
Po = 0.43 GeV cm~ 3 from Salucci et al. (2010), which 
has, however, a large uncertainty, with a typical associated 
error bar of ±0.1 GeV cnT 3 and a possible spread up to 
0.2-0. 7 GeV cm” 3 (Salucci et al. 2010; Cirelli et al. 2011). 
Whether the limits will worsen, or improve, thus awaits a better 
determination of po. To show the effect of the po uncertainty 
on the limits we plot, for illustration, in the top left panel 
of Figures 4 and 5 the uncertainty band (red dotted lines) in 
the 3 a no-background limits which would result from varying 
the local DM density po in the range 0.2-0. 7 GeV cm -3 (we 
conservatively take here the larger scatter to show the maximal 
impact of the uncertainty of po). A similar band, not shown in 
the plot for clarity, would be present for the limits including 
a model of the astrophysical background. The band is likely a 
generous estimate of the uncertainty since the variation of po is 
typically correlated with other properties of the DM halo, such 
as the density profile and the distance of the solar system from 
the Galactic center R$. The uncertainty band shown should just 
be considered an illustration, while a detailed study would be 
required to address the actual uncertainty, which is beyond the 
scope of this work. 

In Figure 4 we also show the regions of the parameter space 
derived in Cirelli et al. (2010) from a DM fit to the Fermi-LAI 
electron/positron data and PAMELA positron fraction data. 
Contours are shown at 95% and 99.999% CL. These regions 
are rescaled down of a factor (0.43/0.3) 2 ~ 2 to take into 
account the different local DM density po used in the two 
works (for the same reason the regions in the decay case are 
rescaled up by a factor of (0.43/0.3) ~ 1.4). We must also 
take into account that different energy losses of local CREs are 
used here when compared to Cirelli et al. (2010). For CREs 
with E ^ 1 TeV diffusion can be approximately neglected so 
that, from the diffusion loss equation, the steady state CRE flux 
scales approximately linearly with the energy loss timescale 
r (see Bergstrom et al. 2009; Borriello et al. 2009a). The 
PAMELA/ Fermi regions, instead, scale as the inverse of the 
local CRE flux. Beyond the different local magnetic field and 
ISRF, a major difference with the analysis in Cirelli et al. (2010) 
is that they neglect Klein-Nishina attenuation effects in the 
IC cross section, which are instead taken into account in our 
GALPR0P calculations. These effects make the energy losses of 
CREs with E ^ 1 TeV almost negligible on the optical part 
of the ISRF, and slightly decrease the energy losses due the 
infrared part. CMB losses are instead unaffected. As a result, 
the energy loss timescale, and thus the CRE flux, increases 
(the PAMELA /Fermi regions will be pushed downward). On 
the other hand, we use a local magnetic field of 5 pG as opposed 
to the 3 pG used in Cirelli et al. (2010), and this has the 
contrary effect, decreasing the energy losses timescale (which 
scales as f?~ 2 ) and thus again pushing up the PA MELA / Fermi 
confidence regions. Overall, to correctly derive the positions 
of the PA M L L A / Fe rm i confidence regions, the uncertainties 
in the local ISRE and magnetic field should be taken into 
account and marginalized away, which we leave for a follow- 
up work. Here, from the above considerations, we estimate 
that the location of the PAMELA/ Fern li regions has a further 
uncertainty of a factor ~2, so that they can possibly touch 
the exclusion limits. Thus, we cannot robustly rule out the 
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Figure 4. Upper limits on the velocity averaged DM annihilation cross section including a model of the astrophysical background compared with the limits obtained 
with no modeling of the background. Upper panel: limits on models in which DM annihilates into bb, for a DM distribution given by the NFW distribution (left) and 
isothermal distribution (right). In the left panel we also add an uncertainty band (red dotted lines) in the 3 cr no-background limits which would result from varying 
the local DM density po in the range 0.2-0. 7 GeV cm -3 . A similar band, not shown in the plot for clarity, would be present for the limits including a model of the 
astrophysical background (see discussion in the text). The horizontal line marks the thermal decoupling cross section expected for a generic WIMP candidate. Middle 
panel: upper limits for DM annihilation to p + p~ . Lower panel: the same, for DM annihilation to t + t . The region excluded by the analysis with no model of the 
astrophysical background is indicated in light blue, while the additional region excluded by the analysis with a modeling of the background is indicated in light green. 
The regions of parameter space which provide a good fit to PAMELA (Adriani et al. 2009a, purple) and Fermi-LXT (Abdo et al. 2009, blue) CR electron and positron 
data are shown, as derived in Cirelli et al. (2010) and are scaled by a factor of 0.5, to account for different assumptions on the local DM density (see the text for more 
details). 

(A color version of this figure is available in the online journal.) 


DM annihilation interpretation of the PAMELA / Fermi CRs 
anomalies, although this interpretation is challenged. Finally, 
we note that the PAMELA region below ~200 GeV is now 
disfavored by the new positron measurements with the LAT 


(Ackermann et al. 2012a), which indicate that the positron 
fraction continues to rise to this energy. 

It should be noted that the above conclusions are not affected 
by the uncertainty in po since both the derived constraints 


13 


The Astrophysical Journal, 761:91 (lBpp), 2012 December 20 


Ackermann et al. 


X bb, NFW 

10 28 = 



— 

jq 25 _ w/o background modeling 

= constrained free source fits 


. , 3f r, p 0 =0.43 GeV cm 3 

i n 24 L 

l 5m, p 0 =0.43 GeV cnr 3 

- 3 o', w/o background modeling, po=0.2-0.7 GeV cm~ 3 

10 23 = 

10 1 0 2 1 0 3 1 0 4 

m [GeVl 


X — > bb, ISO 



* — > pV. NFW 

10 28 [ 



10 23 g 

10 10 2 10 3 10 4 
m [GeV] 


* -> pV, ISO 


10 2! 



10 : 


10 


10 2 


10 3 


I0 + 


m [GeV] 


NFW 


. rV, ISO 



10 - 


10 


10 2 


10 3 


10 4 


m [GeV] 



10 - 


10 


10 2 


10 3 


10 4 


m [GeV] 


Figure 5. Lower limits on the lifetime of decaying DM. The panel structure is the same as in Figure 4. In the top left panel we also add an uncertainty band (red dotted 
lines) in the 3<r no-background limits which would result from varying the local DM density po in the range 0.2-0. 7 GeV cm -3 . A similar band, not shown in the 
plot for clarity, would be present for the limits including a model of the astrophysical background (see discussion in the text). The regions of parameter space which 
provide a good fit to PAMELA (Adriani et al. 2009a, purple) and Fermi - LAT Abdo et al. (2009, blue) CR electron and positron data are shown as derived in (Cirelli 
et al. 2010) and are scaled by a factor of 1.4, to account for different assumptions on the local DM density (see the text for more details). 

(A color version of this figure is available in the online journal.) 


and the region of parameter space compatible with the DM 
interpretation of the CR anomalies scale in the same way with 
Po- The same is true also for the constrains on decaying DM. 

Constraints for the case of decaying DM are shown in 
Figure 5. The interpretation of the PAMELA/ Fermi CR features 
in terms of decaying DM is not ruled out in this analysis. The 


limits are stronger than the ones derived in similar analyses 
performed without background modeling (Papucci & Strumia 
2010; Cirelli et al. 2010) and slightly improved over the ones 
derived from observation of Galaxy clusters (Huang et al. 2012). 
They are comparable to the limits derived from the comparison 
with the IGB (Cirelli et al. 2010; Huang et al. 2012). 
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Figure 6. Best-fit e and p CRSDs and errors, marginalized over the remaining parameters and over the various DM models considered. The pulsar distribution from 
Yusifov & Kii§uk (2004) is also shown for comparison. The source distribution is zeroed within 3 kpc of the Galactic center; see the text. 

(A color version of this figure is available in the online journal.) 


In Figure 6 we show the CR source distributions derived in 
the fit (the coefficients c?’ p ) and the uncertainty of them. Each 
bin is treated as a parameter for which the profile likelihood 
is built by marginalizing over all the other parameters (linear 
and nonlinear) and DM models and 3er and 5er uncertainties are 
derived in the usual way. Again, since the Galactic plane and 
the full energy range is not included in the fit, in interpreting the 
figure the same caveats as for Figure 3 apply. For example, the 
fitted /?CRSD increases in the outer Galaxy: although protons 
from the outer Galaxy can propagate to the inner Galaxy 
(while few also propagate from the other side of the Galactic 
center) and we thus have some sensitivity to them in our ROI, 
clearly, a proper statement on this feature would require to fit 
the outer Galaxy itself. Further caveats in interpreting these 
results include the fact that the CRSDs also depend on the 
ISM distribution and ISRF distribution, which are themselves 
affected by large uncertainties. Finally, the CRSDs also depend 
on the chosen propagation setup. We do not attempt here to 
assess the systematic uncertainties related to these issues (see 
the next section, however, for some discussion regarding the 
ISRF uncertainties). 

Overall, however, it can be noted that the eCRSD and /?CRSD 
are in reasonable agreement with each other and with typical 
SNR or pulsars models, except for the rise in the outer Galaxy 
which is not predicted in these models. The plots are very similar 
when no DM contribution is assumed, which is expected since 
DM gives a subdominant contribution to the gamma-ray signal 
3 kpc away from the GC. 

An interesting point to comment on is how the DM constraints 
are affected if some of the nuisance parameters, in particular 
the CRSDs, are fixed to benchmark choices taken from the 
literature, instead of being marginalized away. We checked the 
resulting DM constraints for cases in which we fix the CRSDs 
(the eCRSD and /?CRSD are taken to be equal) to three different 
common choices, namely SNR (Case & Bhattacharya 1998), 
pulsar (Yusifov & Kiigiik 2004), and a simple Gaussian model 
(centered at the Galactic center, with a width of 4.5 kpc) while 
the rest of the parameters are marginalized away in the same way 
as in the main analysis. We find that the SNR case gives slightly 
worse limits (20%-30%) while pulsar limits are very similar to 
the results of the main analysis (in agreement with the fact that 


the our best-fit CRSDs are close to the pulsar distribution), and 
Gaussian limits are a factor of ^2 better than the pulsar case. 
This last case is understood in light of the fact that a Gaussian 
CRSD, being peaked at the Galactic center, forces the inner 
Galaxy gamma-ray signal to be explained entirely in terms of 
ordinary CR sources, leaving little gamma-ray flux to DM and 
thus giving better limits. This is interesting and might indicate 
that the DM constraints can become better if independent robust 
constraints on the CRSD become available. On the other hand, 
as the results of the main analysis show, the fitted CRSD does 
not favor a Gaussian distribution (at least at a qualitative level; 
we have not performed a quantitative comparison for the reasons 
explained above). For the time being, thus, the use of CRSD as a 
nuisance parameter seems the best approach, which, at the same 
time, leaves freedom in the fit to explore the degeneracy with 
DM and limit the fit to explore CRSDs which are in reasonable 
agreement with the (gamma-ray) data. 

Finally, we show in Figure 7 the counts map in our ROI, 
together with the model prediction for a model close to the 
best fit ( Zh = 10, kpc Ye , 2 = 2.3, and d2HI = 0.0140 x 
10- 20 mag cm 2 ) and its residuals when DM (of mass m x = 
150 GeV and annihilating into bb) is included or not in the fit. 
It can be seen that the residuals are mostly flat, meaning that the 
model (and the models close to the minimum) is a reasonable 
fit to the data. A few features are however present, like the 
excess in the vicinity of (/, b ) ~ (—45, 10) and ~(7, —15) 
which seem to be related to the low-latitude tip of Loop I and 
to the low-latitude part of the south lobe/bubble, respectively. 
The two prominent negative residuals near (Z, b) ~ (—15, 5) and 
~(20, —10), instead, approximately contour the lobes and thus 
seem to be an artifact of the fit to compensate for this missing 
component. Gas misplaced in incorrect annuli also could be an 
alternative explanation. 

We also show the point-source mask used based on the 1FGL 
catalog and, for comparison, the mask based on the 2FGL 
catalog (Nolan et al. 2012) and the residuals using this mask. 
Overall, it can be seen that the 2FGL mask covers few point 
sources which are apparent in the residuals with the 1FGL mask. 
The large-scale features in the residuals are however unchanged, 
apart from a small part of Loop I near (Z, b) — (—45, 10), which 
is resolved into sources. 
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Figure 7. Counts map of the ROI that we consider (upper left panel), model prediction for a model (without DM) close to the best-fit (zh = 10 kpc, = 2.3, and 
d2HI = 0.0140 x 10 20 mag cm 2 ) parameter region (upper right), and residuals in units of a for the same model (second row left) and when DM (of mass m x = 
150 GeV and annihilating into bb) is also included in the fit (second row right). Third row: same as second row but with 2FGL point sources masked instead of 1FGL. 
Fourth row: 1FGL mask (left) and 2FGL mask (right). The model and data counts and the residuals have been smoothed with a 1 .25 Gaussian filter. The point sources 
mask in the residuals have been applied before and after the smoothing. 

(A color version of this figure is available in the online journal.) 


9. DISCUSSION ON MODEL UNCERTAINTIES 

In deriving our limits above we have taken into account 
many possible uncertainties like the ones in the e, /tCRSDs, 
in Zh i the electron index, and the dust to gas ratio. We check 
below the importance of further uncertainties which we have 
not considered explicitly in our scan. 

An important component for which there is still a consider- 
able uncertainty is the ISRF. In particular, the ISRF in the inner 
Galaxy is quite uncertain and the default model we used could 
be a substantial underestimate of the true one in this region. Very 
different ISRFs would affect the propagation of CREs through 
energy losses and this could be especially relevant for the DM 
models in which the IC component is important and provide 
strong constraints, like p channels. Modes dominated by prompt 
radiation, like b and r should, instead, not be significantly sub- 
ject to uncertainties in the ISRF. To perform an explicit check 
we repeated our entire analysis using a different ISRF model 
(T. A. Porter 2012, private communication), which has the bulge 
component increased by a factor of 10 (see Porter et al. 2008 
for a detailed definition), which implies an overall increase in 
the inner Galaxy of a factor of two. The DM limits with this 
enhanced ISRF were, however, not appreciably affected. We 
verified that the enhanced ISRF produces an enhanced IC com- 
ponent, but only within a few degrees of the Galactic center, 


thus not affecting the fit in our ROI. It also should be stressed 
that a more intense ISRF implies more IC emission for the DM 
IC too, so that assuming a lower ISRF gives conservative limits. 
Finally, an ISRF lower than the one assumed here is also pos- 
sible, as the results obtained in Ackermann et al. (2012b, see 
their Figure 1 1) for the CRSD following the pulsar distribution 
seems to indicate. However, the “ISRF normalization” reported 
in Ackermann et al. (2012b) is more precisely a proxy for a com- 
bination of ISRF intensity, normalization of the CRE spectrum, 
and halo size, so that alternative explanations are possible. 

We also checked more systematically other sources of un- 
certainties, but in a more simplified setup: we set a particular 
model as reference and then we varied each parameter one at a 
time, keeping the others fixed, and for each case we calculate 
the percentage variation in DM limits for selected DM models. 
We vary the parameters derived from the CR fit, va, y P ,i, Y P , 2, 
Pbr.p, and the (Do, Zh) relation, within the uncertainty ranges 
derived in Ackermann et al. (20 1 2b) enlarging it by a factor of 
~2 to take into account possible systematic uncertainties (the 
errors quoted in Ackermann et al. 2012b are statistical only). We 
also include in the list of the tested parameters the ones which 
are included in our model scan (CRSD, d2HI, y e 2 , (Do, z/,)) to 
allow for a direct comparison. The following set of parameters, 
which lie close to the best fit of our analysis, was chosen for 
the reference model: Va =36 km s -1 , Do = 5.010 28 cm 2 s _1 , 
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Table 3 

Relative Variation |<5a/a|[%] of the Limits on the DM Velocity Averaged Annihilation Cross-section Derived in This Work 
with Respect to Changes in the Underlying Astrophysical Diffuse Emission Model 


Parameter 

\&a/a\ [%], bb 

\&o/a\ [%], mV" 

v a [ 30; 36; 45] km s _1 

[6; 0; 11] 

[ 4.; 0; 9] 

Y P , i [1.8; 1.9; 2;] 

[ 1.0; 0; 2.5] 

[1.5; 0; 2.0] 

Y P , 2 [ 2.35; 2.39; 2.45] 

[ 2.5; 0; 1.5] 

[2.5; 0; 1.5] 

p br , p [ 10; 11.5; 12.5] GV 

[ 0.5; 0; 1.0] 

[0.9; 0; 1.5] 

d2HI [ 0.0110, 0.0140; 0.0170] 10~ 20 mag cm 2 

[3; 0; 12] 

[ 3;0; 9] 

Ye, 2 [ 2.0; 2.45; 2.6] 

[ 17; 0; 7] 

[ 18; 0; 5] 

(D 0 , z h ) [ (5.0e28, 4); (7.1e28, 10)] cmV 1 

[0; 10] 

[0; 7] 

CRSD [ SNR; Pulsar] 

[ 0; 61] 

[ 0; 59] 

KRA(S = 0.5); KOL(<5 = 0.3); PD(S = 0.6) 

[ 4.0; 0; 3.0] 

[1.0; 0; 5] 

V c [0; 20] km s” 1 

[ 0; 6] 

[0; 4] 

GMF [ Conf 1, Conf 2] 

[0; 3] 

[0; 8] 


Notes. The table shows the relative variation for selected DM models ( bb and channel, for a 150 GeV DM) in 
a simplified setup when only one parameter is varied at a time. Each row corresponds to the indicated parameter. The 
bold values correspond to the reference value. 


Zh = 4 kpc, 8 = 0.3, Yp, 1 = 1-9, y p , 2 = 2.39, p br , p = 11.5, 
Ye , 2 = 2.45, d2HI = 0.014 xlO" 20 mag cm 2 , CRSD = SNR, 
and V c = 0 km s" 1 . Results are shown in Table 3. 

We can see that CR parameters such as Va and Yp, t, Yp. 2, 
Pbr.p , V c , and even different gas maps have very low (<10%) 
impact on the DM limits. The table confirms that Ye ,2 and the 
CRSD (which we fix here to be the same for protons and CREs) 
are the main parameters degenerate with DM and thus affecting 
the limits the most (up to 60%). The diffusion constant Dq is 
tightly correlated to the halo height Zh ■ Therefore we vary the 
parameter pair (Dq, Zh) instead of the single parameters, using 
their relation derived from the fit to the CR data described in 
Section 2. Nonetheless, the combination of Dq and Zh is included 
individually in the parameter scans of the previous section used 
for the main results. As an additional check of the effect of the 
CR propagation parameters on the DM limits, we find DM limits 
in three theoretical CR propagation setups: plain diffusion (PD, 
characterized by index of diffusion of 8 = 0.6), Kraichanian 
(KRA, 8 = 0.5), and Kolmogorov (KOL, 8 = 0.3). In these 
cases, the rest of the CR propagation parameters are found from 
the best fit to the CR data following the method described in 
Ackermann et al. (2012b) and di Bernardo et al. (2011). These 
fits to CR are performed again without any DM component. We 
find that DM limits in these three CR diffusion setups are also 
barely affected, in particular when compared to the effect of the 
CR source distribution, as shown in Table 3. 

Finally, we also consider an alternative configuration of the 
GMF. The reference one is the default configuration used in 
GALPR0P with an exponential profile in R and z and length 
scales of 10 kpc in R and 2 kpc in z, normalized to 5 /rG locally 
(Conf 1). The alternative configuration we tested has in addition 
a further component of constant 100 //G intensity within 0.4 kpc 
from the Galactic center, as motivated by a recent work 
Crocker et al. (2010; Conf 2). This alternative configuration 
also produces changes in the limits of less than 10%. 

10. SUMMARY 

In this work we constrain the contribution to the diffuse 
gamma-ray emission from DM annihilating or decaying in the 
Milky Way halo, based on the Fenni-LAT gamma-ray data. We 
first present the most conservative limits on DM assuming that 
all FAT photons from the halo are produced by annihilating/ 
decaying WIMPs. Then, based on our current best knowledge 


of the Galactic diffuse emission (Ackermann et al. 2012b), we 
use GALPR0P to model the astrophysical diffuse background, 
and, using a profile likelihood approach, we explore the effects 
of various poorly constrained parameters in the modeling of 
the astrophysical background, e.g., the diffusive halo height, 
the CR source distribution, the index of the electron injection 
spectrum, and the dust to gas ratio in order to get more robust 
constraints. We also remove astrophysical CR sources within 
3 kpc from the Galactic center so that any potential astrophysical 
contribution in this region is attributed to DM, resulting in 
more conservative constraints. Overall, rather than being due 
to residual astrophysical model uncertainties, the remaining 
major uncertainties in the DM constraints from the halo region 
come from the modeling of the DM signal itself. The main 
uncertainty is in the normalization of the DM profile, which is 
fixed through the local value of the DM density. We use the 
recent determination pq = 0.43 GeV cm -3 from Salucci et al. 
(2010), which has, however, a large uncertainty, with values in 
the range 0. 2-0.7 GeV cm" 3 still viable. A large uncertainty in 
Pq is particular important for annihilation constraints since they 
scale like pj, while for constraints on decaying DM the scaling is 
only linear. A less important role is played by the uncertainties 
in the DM profile, since in the halo region different profiles 
predict similar DM densities. When using the lowest allowed 
DM density used in literature of po = 0.2 GeV cm" 3 (see, e.g., 
Salucci et al. 2010), our limits worsen by a factor 4 (2) for 
annihilating (decaying) DM, as shown in Figures 4 and 5. A 
better determination of the local DM density, as well as of the 
parameters determining the global structure of the DM halo, is 
therefore of the utmost importance for reducing the uncertainties 
related to DM constraints from DM halo, but it is beyond the 
scope of this paper and is the subject of dedicated studies. 

Bearing this in mind, the limits we obtain are competitive with 
complementary probes of DM like dwarfs, clusters of galaxies, 
or recent constraints obtained from CMB observations for DM 
models with prompt spectra, and significantly improve over 
these studies for DM models with significant IC contribution 
such as DM annihilating into p + /i". The limits we derive for 
leptonic models challenge the interpretation of the PAMEFA 
and Fermi CR anomalies as annihilation of DM in the Galactic 
halo, while they are not constraining enough to exclude the 
interpretation in terms of decaying DM. We note that this last 
conclusion is not affected by the uncertainty in po since both the 
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derived constraints and the region of parameter space compatible 
with the DM interpretation of the CR anomalies scale in the same 
way with po- 

An obvious improvement of this analysis would be a full scan 
of the CR parameters from a simultaneous fit to both gamma 
and CR data. An effort in this direction is currently ongoing and 
will be reported in future publications. 
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